Preparations

Load the necessary libraries

library(rstanarm)   #for fitting models in STAN
library(brms)       #for fitting models in STAN
library(coda)       #for diagnostics
library(bayesplot)  #for diagnostics
library(ggmcmc)     #for MCMC diagnostics
library(DHARMa)     #for residual diagnostics
library(rstan)      #for interfacing with STAN
library(emmeans)    #for marginal means etc
library(broom)      #for tidying outputs
library(tidybayes)  #for more tidying outputs
library(ggeffects)  #for partial plots
library(tidyverse)  #for data wrangling etc
library(broom.mixed)#for summarising models
library(ggeffects)  #for partial effects plots
theme_set(theme_grey()) #put the default ggplot theme back

Scenario

Polis et al. (1998) were intested in modelling the presence/absence of lizards (Uta sp.) against the perimeter to area ratio of 19 islands in the Gulf of California.

Uta lizard

Format of polis.csv data file

ISLAND	RATIO	PA
Bota	15.41	1
Cabeza	5.63	1
Cerraja	25.92	1
Coronadito	15.17	0
..	..	..

ISLAND	Categorical listing of the name of the 19 islands used - variable not used in analysis.
RATIO	Ratio of perimeter to area of the island.
PA	Presence (1) or absence (0) of Uta lizards on island.

The aim of the analysis is to investigate the relationship between island parimeter to area ratio and the presence/absence of Uta lizards.

Read in the data

polis = read_csv('../public/data/polis.csv', trim_ws=TRUE)
glimpse(polis)

## Rows: 19
## Columns: 3
## $ ISLAND <chr> "Bota", "Cabeza", "Cerraja", "Coronadito", "Flecha", "Gemelose…
## $ RATIO  <dbl> 15.41, 5.63, 25.92, 15.17, 13.04, 18.85, 30.95, 22.87, 12.01, …
## $ PA     <dbl> 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1

head(polis)

str(polis)

## tibble [19 × 3] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ ISLAND: chr [1:19] "Bota" "Cabeza" "Cerraja" "Coronadito" ...
##  $ RATIO : num [1:19] 15.41 5.63 25.92 15.17 13.04 ...
##  $ PA    : num [1:19] 1 1 1 0 1 0 0 0 0 1 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   ISLAND = col_character(),
##   ..   RATIO = col_double(),
##   ..   PA = col_double()
##   .. )

Exploratory data analysis

The individual responses (\(y_i\), observed presence/absence of Uta lizards) are each expected to have been independently drawn from Bernoulli (or binomial) distributions (\(\mathcal{Bin}\)). These distributions represent all the possible presence/absences we could have obtained at the specific (\(i^th\)) level of island perimeter to area ratio. Hence the \(i^th\) presence/absence observation is expected to have been drawn from a binomial distribution with a probability of \(\mu_i\) and size of (\(n=1\)).

The expected probabilities are related to the linear predictor (intercept plus slope associated with perimeter to area ratio) via a logit link.

We need to supply priors for each of the parameters to be estimated (\(\beta_0\) and \(\beta_1\)). Whilst we want these priors to be sufficiently vague as to not influence the outcomes of the analysis (and thus be equivalent to the frequentist analysis), we do not want the priors to be so vague (wide) that they permit the MCMC sampler to drift off into parameter space that is both ilogical as well as numerically awkward.

As a starting point, lets assign the following priors:

\(\beta_0\): Normal prior centered at 0 with a variance of 10
\(\beta_1\): Normal prior centered at 0 with a variance of 1

Note, when fitting models through either rstanarm or brms, the priors assume that the predictor(s) have been centered and are to be applied on the link scale. In this case the link scale is an identity.

Model formula: \[ \begin{align} y_i &\sim{} \mathcal{Bin}(n, p_i)\\ ln\left(\frac{p_i}{1-p_i}\right) &= \beta_0 + \beta_1 x_i\\ \beta_0 &\sim{} \mathcal{N}(0,10)\\ \beta_1 &\sim{} \mathcal{N}(0,1)\\ \end{align} \]

Fit the model

frequentist

summary(glm(PA ~ RATIO, data=polis, family=binomial()))

## 
## Call:
## glm(formula = PA ~ RATIO, family = binomial(), data = polis)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6067  -0.6382   0.2368   0.4332   2.0986  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   3.6061     1.6953   2.127   0.0334 *
## RATIO        -0.2196     0.1005  -2.184   0.0289 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 26.287  on 18  degrees of freedom
## Residual deviance: 14.221  on 17  degrees of freedom
## AIC: 18.221
## 
## Number of Fisher Scoring iterations: 6

rstanarm

Using default priors

In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.

polis.rstanarm = stan_glm(PA ~ RATIO, data=polis,
                          family=binomial(), 
                         iter = 5000, warmup = 1000,
                         chains = 3, thin = 5, refresh = 0)

In the above:

the formula and data arguments should be familiar as they are the same as for many models in R
iter: indicates the number of MCMC iterations to perform per chain
warmup: indicates how many of the initial MCMC samples to exclude. During the warmup stage, the MCMC sampler is tuning its sampling operations to allow it to determine the best way to create the chain. While it is in this discovery phase, the estimates it produces might be unrepresentative.
chains: indicates the number of independent chains to run. The more chains, the more random starts. The point of multiple random starts is to mimimise the potential for the starting point to dictate the outcome. It is possible that a chain might get stuck sampling a single feature of the posterior and not discover additional features. Additional random starting points should hopefully elivate this. Either 3 or 4 is typical.
thin: indicates the rate at which the MCMC samples should be thinned in order to eliminate autocorrelation between MCMC samples. When the step lengths between MCMC samples are small, the resulting parameter estimates will be autocorrelated (since samples close by in the chain will be more similar than those far appart). The rate of thinning necessary will depend on the degree of autocorrelation. For stan models, a thinning rate of 5 is a good starting point.
refresh: indicates how often the terminal display should be updated. When running a large model, this can be useful as it provides a form of progress. However, when the routine is within an Rmardown document, it just results in the output including each line of progress and a lot of space is taken up unnecessarily

Having allowed rstanarm to formulate its own weakly informative priors, it is a good idea to explore what they are. Firstly, out of curiosity, it might be interesting to see what it has chosen. However, more importantly, we need to be able to document what the priors were and the rstanarm development team make it very clear that there is no guarantee that the default priors will remain the same into the future.

prior_summary(polis.rstanarm)

## Priors for model 'polis.rstanarm' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = 0, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 0, scale = 0.14)
## ------
## See help('prior_summary.stanreg') for more details

This tells us:

for the intercept, it is using a normal prior with a mean of 0.53 and a standard deviation of 2.5. This mean value is based on the mean of the response variable (rounded up). The 2.5 is used for all intercepts. The scale of the prior is then adjusted proportional to the variability in the response, by multiplying 2.5 by the standard deviation of the response (and rounding up).

mean(polis$PA)

## [1] 0.5263158

sd(polis$PA)

## [1] 0.5129892

2.5*sd(polis$PA)

## [1] 1.282473

for the coefficients (in this case, just the slope), the default prior is a normal prior centered around 0 with a standard deviation of 2.5. For binomial models, this is then adjusted for the scale of the data by dividing the 2.5 by the standard deviation of the predictor (then rounded).

2.5/sd(polis$RATIO)

## [1] 0.1429651

there is no auxillary parameter and thus there is no auxillary prior ### Assessing priors

One way to assess the priors is to have the MCMC sampler sample purely from the prior predictive distribution without conditioning on the observed data. Doing so provides a glimpse at the range of predictions possible under the priors. On the one hand, wide ranging preditions would ensure that the priors are unlikely to influence the actual preditions once they are conditioned on the data. On the other hand, if they are too wide, the sampler is being permitted to traverse into regions of parameter space that are not logically possible in the context of the actual underlying ecological context. Not only could this mean that illogical parameter estimates are possible, when the sampler is traversing regions of paramter space that are not supported by the actual data, the sampler can become unstable and have difficulty.

We can draw from the prior predictive distribution instead of conditioning on the response, by updating the model and indicating prior_PD=TRUE. After refittin the model in this way, we can plot the predictions to gain insights into the range of predictions supported by the priors alone.

polis.rstanarm1 <- update(polis.rstanarm,  prior_PD=TRUE)

ggemmeans(polis.rstanarm1,  ~RATIO) %>% plot(add.data=TRUE)

Conclusions:

we see that the range of predictions is very wide and the slope could range from strongly negative to strongly positive.

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):

\(\beta_0\): normal centered at 0.5 with a standard deviation of 10
\(\beta_1\): normal centered at 0 with a standard deviation of 1

I will also overlay the raw data for comparison.

polis.rstanarm2= stan_glm(PA ~ RATIO, data=polis,
                          family=binomial(), 
                          prior_intercept = normal(0.5, 2, autoscale=FALSE),
                          prior = normal(0, 0.2, autoscale=FALSE),
                          prior_PD=TRUE, 
                          iter = 5000, warmup = 1000,
                          chains = 3, thin = 5, refresh = 0
                          )

ggemmeans(polis.rstanarm2,  ~RATIO) %>%
  plot(add.data=TRUE)

Now lets refit, conditioning on the data.

polis.rstanarm3= update(polis.rstanarm2,  prior_PD=FALSE)

Plotting prior and posterior

posterior_vs_prior(polis.rstanarm3, color_by='vs', group_by=TRUE,
                   facet_args=list(scales='free_y'))

Conclusions:

in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.

ggemmeans(polis.rstanarm3,  ~RATIO) %>% plot(add.data=TRUE)

brms

Using default priors

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

polis.brm = brm(bf(PA|trials(1) ~ RATIO, family=binomial()),
                data=polis,
                iter = 5000, warmup = 1000,
                chains = 3, thin = 5, refresh = 0)

In the above:

the formula and data arguments should be familiar as they are the same as for many models in R. That said, the formual specification of brms can define models that are not possible by most other routines. To facilitate this enhanced functionality, we usually define a brms formula within its own bf() function along with the family (in this case, it is gaussian, which is the default and therefore can be omitted.)
iter: indicates the number of MCMC iterations to perform per chain
warmup: indicates how many of the initial MCMC samples to exclude. During the warmup stage, the MCMC sampler is tuning its sampling operations to allow it to determine the best way to create the chain. While it is in this discovery phase, the estimates it produces might be unrepresentative.
chains: indicates the number of independent chains to run. The more chains, the more random starts. The point of multiple random starts is to mimimise the potential for the starting point to dictate the outcome. It is possible that a chain might get stuck sampling a single feature of the posterior and not discover additional features. Additional random starting points should hopefully elivate this. Either 3 or 4 is typical.
thin: indicates the rate at which the MCMC samples should be thinned in order to eliminate autocorrelation between MCMC samples. When the step lengths between MCMC samples are small, the resulting parameter estimates will be autocorrelated (since samples close by in the chain will be more similar than those far appart). The rate of thinning necessary will depend on the degree of autocorrelation. For stan models, a thinning rate of 5 is a good starting point.
refresh: indicates how often the terminal display should be updated. When running a large model, this can be useful as it provides a form of progress. However, when the routine is within an Rmardown document, it just results in the output including each line of progress and a lot of space is taken up unnecessarily

prior_summary(polis.brm)

##                 prior     class  coef group resp dpar nlpar bound       source
##                (flat)         b                                        default
##                (flat)         b RATIO                             (vectorized)
##  student_t(3, 0, 2.5) Intercept                                        default

This tells us:

for the intercept, it is using a student t (flatter normal) prior with a mean of 0 and a standard deviation of 2.5. These are the defaults used for a binomial model.
for the beta coefficients (in this case, just the slope), the default prior is a inproper flat prior. A flat prior essentially means that any value between negative infinity and positive infinity are equally likely. Whilst this might seem reckless, in practice, it seems to work reasonably well for non-intercept beta parameters.
there is no sigma parameter in a binomial model and therefore there are no additional priors.

Assessing priors

In brms, we can inform the sampler to draw from the prior predictive distribution instead of conditioning on the response, by running the model with the sample_prior='only' argument. Unfortunately, this cannot be applied when there are flat priors (since the posteriors will necessarily extend to negative and positive infinity). Therefore, in order to use this useful routine, we need to make sure that we have defined a proper prior for all parameters.

polis.brm1 = brm(bf(PA|trials(1) ~ RATIO, family=binomial()),
                 data=polis,
                prior=c(
                  prior(normal(0, 10), class='Intercept'), 
                  prior(normal(0, 1), class='b')), 
                sample_prior = 'only', 
                iter = 5000, warmup = 1000,
                chains = 3, thin = 5, refresh = 0)

ggemmeans(polis.brm1,  ~RATIO) %>% plot(add.data=TRUE)

conditional_effects(polis.brm1) %>%  plot(points=TRUE)

Conclusions:

we see that the range of predictions is faily wide and the slope could range from strongly negative to strongly positive.

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):

\(\beta_0\): normal centered at 0 with a standard deviation of 10
\(\beta_1\): normal centered at 0 with a standard deviation of 1

I will also overlay the raw data for comparison.

polis.brm2 = brm(bf(PA|trials(1) ~ RATIO, family=binomial()),
                 data=polis,
                 prior=c(prior(normal(0, 10),  class='Intercept'),
                            prior(normal(0, 1), class='b')
                            ), 
                 sample_prior = 'only', 
                 iter = 5000, warmup = 1000,
                 chains = 3, thin = 5, refresh = 0)

ggemmeans(polis.brm2,  ~RATIO) %>%
  plot(add.data=TRUE)

polis.brm3 <- update(polis.brm2,  sample_prior=TRUE, refresh=0)

Plotting prior and posterior

polis.brm3 %>% get_variables()

##  [1] "b_Intercept"     "b_RATIO"         "prior_Intercept" "prior_b"        
##  [5] "lp__"            "accept_stat__"   "stepsize__"      "treedepth__"    
##  [9] "n_leapfrog__"    "divergent__"     "energy__"

polis.brm3 %>%
  posterior_samples %>%
  select(-`lp__`) %>%
  gather %>%
  mutate(Type=ifelse(str_detect(key, 'prior'), 'Prior', 'b'),
         Class=ifelse(str_detect(key, 'Intercept'),  'Intercept',
               ifelse(str_detect(key, 'sigma'),  'Sigma',  'b'))) %>%
  ggplot(aes(x=Type,  y=value)) +
  stat_pointinterval()+
  facet_wrap(~Class,  scales='free')

Exploring the stan code

standata(polis.brm3)

## $N
## [1] 19
## 
## $Y
##  [1] 1 1 1 0 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1
## 
## $trials
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## $K
## [1] 2
## 
## $X
##    Intercept RATIO
## 1          1 15.41
## 2          1  5.63
## 3          1 25.92
## 4          1 15.17
## 5          1 13.04
## 6          1 18.85
## 7          1 30.95
## 8          1 22.87
## 9          1 12.01
## 10         1 11.60
## 11         1  6.09
## 12         1  2.28
## 13         1  4.05
## 14         1 59.94
## 15         1 63.16
## 16         1 22.76
## 17         1 23.54
## 18         1  0.21
## 19         1  2.55
## attr(,"assign")
## [1] 0 1
## 
## $prior_only
## [1] 0
## 
## attr(,"class")
## [1] "standata"

stancode(polis.brm3)

## // generated with brms 2.14.0
## functions {
## }
## data {
##   int<lower=1> N;  // total number of observations
##   int Y[N];  // response variable
##   int trials[N];  // number of trials
##   int<lower=1> K;  // number of population-level effects
##   matrix[N, K] X;  // population-level design matrix
##   int prior_only;  // should the likelihood be ignored?
## }
## transformed data {
##   int Kc = K - 1;
##   matrix[N, Kc] Xc;  // centered version of X without an intercept
##   vector[Kc] means_X;  // column means of X before centering
##   for (i in 2:K) {
##     means_X[i - 1] = mean(X[, i]);
##     Xc[, i - 1] = X[, i] - means_X[i - 1];
##   }
## }
## parameters {
##   vector[Kc] b;  // population-level effects
##   real Intercept;  // temporary intercept for centered predictors
## }
## transformed parameters {
## }
## model {
##   // likelihood including all constants
##   if (!prior_only) {
##     // initialize linear predictor term
##     vector[N] mu = Intercept + Xc * b;
##     target += binomial_logit_lpmf(Y | trials, mu);
##   }
##   // priors including all constants
##   target += normal_lpdf(b | 0, 1);
##   target += normal_lpdf(Intercept | 0, 10);
## }
## generated quantities {
##   // actual population-level intercept
##   real b_Intercept = Intercept - dot_product(means_X, b);
##   // additionally draw samples from priors
##   real prior_b = normal_rng(0,1);
##   real prior_Intercept = normal_rng(0,10);
## }

MCMC sampling diagnostics

In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.

There are a wide variety of tests that range from the big picture, overal chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.

rstanarm

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

available_mcmc()

## bayesplot MCMC module:
##   mcmc_acf
##   mcmc_acf_bar
##   mcmc_areas
##   mcmc_areas_data
##   mcmc_areas_ridges
##   mcmc_areas_ridges_data
##   mcmc_combo
##   mcmc_dens
##   mcmc_dens_chains
##   mcmc_dens_chains_data
##   mcmc_dens_overlay
##   mcmc_hex
##   mcmc_hist
##   mcmc_hist_by_chain
##   mcmc_intervals
##   mcmc_intervals_data
##   mcmc_neff
##   mcmc_neff_data
##   mcmc_neff_hist
##   mcmc_nuts_acceptance
##   mcmc_nuts_divergence
##   mcmc_nuts_energy
##   mcmc_nuts_stepsize
##   mcmc_nuts_treedepth
##   mcmc_pairs
##   mcmc_parcoord
##   mcmc_parcoord_data
##   mcmc_rank_hist
##   mcmc_rank_overlay
##   mcmc_recover_hist
##   mcmc_recover_intervals
##   mcmc_recover_scatter
##   mcmc_rhat
##   mcmc_rhat_data
##   mcmc_rhat_hist
##   mcmc_scatter
##   mcmc_trace
##   mcmc_trace_data
##   mcmc_trace_highlight
##   mcmc_violin

Of these, we will focus on:

mcmc_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

plot(polis.rstanarm3, plotfun='mcmc_trace')

The chains appear well mixed and very similar

acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

plot(polis.rstanarm3, 'acf_bar')

There is no evidence of autocorrelation in the MCMC samples

Rhat: Rhat is a measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

plot(polis.rstanarm3, 'rhat_hist')

All Rhat values are below 1.05, suggesting the chains have converged.

neff (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

plot(polis.rstanarm3, 'neff_hist')

Ratios all very high.

More diagnostics

plot(polis.rstanarm3, 'combo')

plot(polis.rstanarm3, 'violin')

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

stan_trace(polis.rstanarm3)

The chains appear well mixed and very similar

stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

stan_ac(polis.rstanarm3)

There is no evidence of autocorrelation in the MCMC samples

stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

stan_rhat(polis.rstanarm3)

All Rhat values are below 1.05, suggesting the chains have converged.

stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

stan_ess(polis.rstanarm3)

Ratios all very high.

stan_dens(polis.rstanarm3, separate_chains = TRUE)

ggmcmc

The ggmean package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

polis.ggs <- ggs(polis.rstanarm3)
ggs_traceplot(polis.ggs)

The chains appear well mixed and very similar

gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

ggs_autocorrelation(polis.ggs)

There is no evidence of autocorrelation in the MCMC samples

stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

ggs_Rhat(polis.ggs)

All Rhat values are below 1.05, suggesting the chains have converged.

stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

ggs_effective(polis.ggs)

Ratios all very high.

More diagnostics

ggs_crosscorrelation(polis.ggs)

ggs_grb(polis.ggs)

brms

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

available_mcmc()

## bayesplot MCMC module:
##   mcmc_acf
##   mcmc_acf_bar
##   mcmc_areas
##   mcmc_areas_data
##   mcmc_areas_ridges
##   mcmc_areas_ridges_data
##   mcmc_combo
##   mcmc_dens
##   mcmc_dens_chains
##   mcmc_dens_chains_data
##   mcmc_dens_overlay
##   mcmc_hex
##   mcmc_hist
##   mcmc_hist_by_chain
##   mcmc_intervals
##   mcmc_intervals_data
##   mcmc_neff
##   mcmc_neff_data
##   mcmc_neff_hist
##   mcmc_nuts_acceptance
##   mcmc_nuts_divergence
##   mcmc_nuts_energy
##   mcmc_nuts_stepsize
##   mcmc_nuts_treedepth
##   mcmc_pairs
##   mcmc_parcoord
##   mcmc_parcoord_data
##   mcmc_rank_hist
##   mcmc_rank_overlay
##   mcmc_recover_hist
##   mcmc_recover_intervals
##   mcmc_recover_scatter
##   mcmc_rhat
##   mcmc_rhat_data
##   mcmc_rhat_hist
##   mcmc_scatter
##   mcmc_trace
##   mcmc_trace_data
##   mcmc_trace_highlight
##   mcmc_violin

Of these, we will focus on:

trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

mcmc_plot(polis.brm3, type='trace')

The chains appear well mixed and very similar

acf_bar (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

mcmc_plot(polis.brm3, type='acf_bar')

There is no evidence of autocorrelation in the MCMC samples

rhat_hist: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

mcmc_plot(polis.brm3, type='rhat_hist')

All Rhat values are below 1.05, suggesting the chains have converged.

neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

mcmc_plot(polis.brm3, type='neff_hist')

Ratios all very high.

More diagnostics

mcmc_plot(polis.brm3, type='combo')

mcmc_plot(polis.brm3, type='violin')

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

stan_trace(polis.brm3$fit)

The chains appear well mixed and very similar

stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

stan_ac(polis.brm3$fit)

There is no evidence of autocorrelation in the MCMC samples

stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

stan_rhat(polis.brm3$fit)

All Rhat values are below 1.05, suggesting the chains have converged.

stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

stan_ess(polis.brm3$fit)

Ratios all very high.

stan_dens(polis.brm3$fit, separate_chains = TRUE)

ggmcmc

The ggmean package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

polis.ggs <- ggs(polis.brm3)
ggs_traceplot(polis.ggs)

The chains appear well mixed and very similar

gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

ggs_autocorrelation(polis.ggs)

There is no evidence of autocorrelation in the MCMC samples

stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

ggs_Rhat(polis.ggs)

All Rhat values are below 1.05, suggesting the chains have converged.

stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

ggs_effective(polis.ggs)

Ratios all very high.

More diagnostics

ggs_crosscorrelation(polis.ggs)

ggs_grb(polis.ggs)

Model validation

rstanarm

pp check

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()

## bayesplot PPC module:
##   ppc_bars
##   ppc_bars_grouped
##   ppc_boxplot
##   ppc_data
##   ppc_dens
##   ppc_dens_overlay
##   ppc_ecdf_overlay
##   ppc_error_binned
##   ppc_error_hist
##   ppc_error_hist_grouped
##   ppc_error_scatter
##   ppc_error_scatter_avg
##   ppc_error_scatter_avg_vs_x
##   ppc_freqpoly
##   ppc_freqpoly_grouped
##   ppc_hist
##   ppc_intervals
##   ppc_intervals_data
##   ppc_intervals_grouped
##   ppc_loo_intervals
##   ppc_loo_pit
##   ppc_loo_pit_overlay
##   ppc_loo_pit_qq
##   ppc_loo_ribbon
##   ppc_ribbon
##   ppc_ribbon_data
##   ppc_ribbon_grouped
##   ppc_rootogram
##   ppc_scatter
##   ppc_scatter_avg
##   ppc_scatter_avg_grouped
##   ppc_stat
##   ppc_stat_2d
##   ppc_stat_freqpoly_grouped
##   ppc_stat_grouped
##   ppc_violin_grouped

dens_overlay: plots the density distribution of the observed data (black line) overlayed ontop of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

pp_check(polis.rstanarm3,  plotfun='dens_overlay')

The model draws appear to be consistent with the observed data.

error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. That said, these are not overly interpretable for models involving binary data.

pp_check(polis.rstanarm3, plotfun='error_scatter_avg')

This is not interpretable.

error_scatter_avg_vs_x: this is similar to a regular residual plot and as such should be interpreted as such. Again, this is not interpretable for binary data.

pp_check(polis.rstanarm3, x=polis$RATIO, plotfun='error_scatter_avg_vs_x')

intervals: plots the observed data overlayed ontop of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

pp_check(polis.rstanarm3, x=polis$RATIO, plotfun='intervals')

The modelled preditions are mostly consitent with the observed data. There is really only one exception.

ribbon: this is just an alternative way of expressing the above plot.

pp_check(polis.rstanarm3, x=polis$RATIO, plotfun='ribbon')

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

#library(shinystan)
#launch_shinystan(polis.rstanarm3)

DHARMa residuals

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components ourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

simulated (predicted) responses associated with each observation.
observed values
fitted (predicted) responses (averaged) associated with each observation

preds <- posterior_predict(polis.rstanarm3,  nsamples=250,  summary=FALSE)
polis.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = polis$PA,
                            fittedPredictedResponse = apply(preds, 2, median),
                            integerResponse = TRUE)
plot(polis.resids)

Conclusions:

the simulated residuals do not suggest any issues with the residuals
there is no evidence of a lack of fit.

brms

pp check

available_ppc()

## bayesplot PPC module:
##   ppc_bars
##   ppc_bars_grouped
##   ppc_boxplot
##   ppc_data
##   ppc_dens
##   ppc_dens_overlay
##   ppc_ecdf_overlay
##   ppc_error_binned
##   ppc_error_hist
##   ppc_error_hist_grouped
##   ppc_error_scatter
##   ppc_error_scatter_avg
##   ppc_error_scatter_avg_vs_x
##   ppc_freqpoly
##   ppc_freqpoly_grouped
##   ppc_hist
##   ppc_intervals
##   ppc_intervals_data
##   ppc_intervals_grouped
##   ppc_loo_intervals
##   ppc_loo_pit
##   ppc_loo_pit_overlay
##   ppc_loo_pit_qq
##   ppc_loo_ribbon
##   ppc_ribbon
##   ppc_ribbon_data
##   ppc_ribbon_grouped
##   ppc_rootogram
##   ppc_scatter
##   ppc_scatter_avg
##   ppc_scatter_avg_grouped
##   ppc_stat
##   ppc_stat_2d
##   ppc_stat_freqpoly_grouped
##   ppc_stat_grouped
##   ppc_violin_grouped

dens_overlay: plots the density distribution of the observed data (black line) overlayed ontop of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

pp_check(polis.brm3,  type='dens_overlay')

The model draws appear to be consistent with the observed data.

error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. Note, this is not really that useful for models that involve a binomial response

pp_check(polis.brm3, type='error_scatter_avg')

This is not really interpretable

error_scatter_avg_vs_x: this is similar to a regular residual plot and as such should be interpreted as such. This is not easy to interpret for binomial responses.

pp_check(polis.brm3, x='RATIO', type='error_scatter_avg_vs_x')

intervals: plots the observed data overlayed ontop of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

pp_check(polis.brm3, x='RATIO', type='intervals')

ribbon: this is just an alternative way of expressing the above plot.

pp_check(polis.brm3, x='RATIO', type='ribbon')

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

#library(shinystan)
#launch_shinystan(polis.brm3)

DHARMa residuals

We need to supply:

simulated (predicted) responses associated with each observation.
observed values
fitted (predicted) responses (averaged) associated with each observation

preds <- posterior_predict(polis.brm3,  nsamples=250,  summary=FALSE)
polis.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = polis$PA,
                            fittedPredictedResponse = apply(preds, 2, median),
                            integerResponse = TRUE)
plot(polis.resids)

Conclusions:

the simulated residuals do not suggest any issues with the residuals
there is no evidence of a lack of fit.

Partial effects plots

rstanarm

ggpredict

polis.rstanarm3 %>% ggpredict() %>% plot(add.data=TRUE)

## $RATIO

ggemmeans

polis.rstanarm3 %>% ggemmeans(~RATIO) %>% plot(add.data=TRUE)

fitted_draws

polis.rstanarm3 %>% fitted_draws(newdata=polis) %>%
 median_hdci() %>%
 ggplot(aes(x=RATIO, y=.value)) +
 geom_ribbon(aes(ymin=.lower, ymax=.upper), fill='blue', alpha=0.3) + 
 geom_line()

brms

conditional_effects

polis.brm3 %>% conditional_effects()

polis.brm3 %>% conditional_effects(spaghetti=TRUE,nsamples=200)

ggpredict

polis.brm3 %>% ggpredict() %>% plot(add.data=TRUE)

## $RATIO

ggemmeans

polis.brm3 %>% ggemmeans(~RATIO) %>% plot(add.data=TRUE)

fitted_draws

polis.brm3 %>% fitted_draws(newdata=polis) %>%
 median_hdci() %>%
 ggplot(aes(x=RATIO, y=.value)) +
 geom_ribbon(aes(ymin=.lower, ymax=.upper), fill='blue', alpha=0.3) + 
 geom_line()

Model investigation

rstanarm

rstanarm captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

summary(polis.rstanarm3)

## 
## Model Info:
##  function:     stan_glm
##  family:       binomial [logit]
##  formula:      PA ~ RATIO
##  algorithm:    sampling
##  sample:       2400 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 19
##  predictors:   2
## 
## Estimates:
##               mean   sd   10%   50%   90%
## (Intercept)  4.6    2.0  2.3   4.4   7.2 
## RATIO       -0.3    0.1 -0.4  -0.3  -0.1 
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 0.5    0.1  0.4   0.5   0.6  
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.0  1.0  2259 
## RATIO         0.0  1.0  2275 
## mean_PPD      0.0  1.0  2212 
## log-posterior 0.0  1.0  2216 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Conclusions:

in the Model Info, we are informed that the total MCMC posterior sample size is 2400 and that there were 19 raw observations.
the estimated intercept (expected presence/absence when perimeter to area ratio is 0) is 4.61. This is the mean of the posterior distribution for this parameter. If we back-transform this to the odds-ratio scale, this becomes 100.71. When the perimeter to area ratio is 0, Uta lizards are 100.71 times more likely to be present than absent. We can also backtransform to the probability scale. When the perimeter to area ratio is 0, the probability of Uta lizards being present is 0.99
the estimated slope (rate at which presence/absence changes per 1 unit change in perimeter to area ratio), is -0.28 (mean) or -0.43 (median) with a standard deviation of 0. The 90% credibility intervals indicate that we are 90% confident that the slope is between -0.28 and -0.27 - e.g. there is a significant positive trend. When backtransformed onto the odds scale, we see that for a one unit increase in perimeter to area ratio, the odds of Uta lizards being present declines by 0.75
Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

tidyMCMC

tidyMCMC(polis.rstanarm3$stanfit, estimate.method='median',  conf.int=TRUE,
         conf.method='HPDinterval',  rhat=TRUE, ess=TRUE)

Conclusions:

the estimated intercept (expected presence/absence when perimeter to area ratio is 0) is 4.41. This is the median of the posterior distribution for this parameter. If we back-transform this to the odds-ratio scale, this becomes 82.53. When the perimeter to area ratio is 0, Uta lizards are 82.53 times more likely to be present than absent. We can also backtransform to the probability scale. When the perimeter to area ratio is 0, the probability of Uta lizards being present is 0.99
the estimated slope (rate at which presence/absence changes per 1 unit change in perimeter to area ratio), is -0.27 (median) with a standard error of 0.12. The 95% credibility intervals indicate that we are 95% confident that the slope is between -0.52 and -0.08 - e.g. there is a significant positive trend. When backtransformed onto the odds scale, we see that for a one unit increase in perimeter to area ratio, the odds of Uta lizards being present declines by 0.77
Rhat and number of effective samples for each parameter are also provided and all look good.

gather_draws

polis.draw <- polis.rstanarm3 %>% gather_draws(`(Intercept)`, RATIO)
## OR via regex
polis.draw <- polis.rstanarm3 %>% gather_draws(`.Intercept.*|RATIO.*`,  regex=TRUE)
polis.draw

We can then summarise this

polis.draw %>% median_hdci

polis.rstanarm3 %>% 
  gather_draws(`(Intercept)`, RATIO) %>% 
  ggplot() + 
  stat_halfeye(aes(x=.value,  y=.variable)) +
  facet_wrap(~.variable, scales='free')

We could alternatively express the parameters on the odds (odds ratio) scale.

polis.rstanarm3 %>% 
  gather_draws(`(Intercept)`, RATIO) %>%
  group_by(.variable) %>%
  mutate(.value=exp(.value)) %>%
  median_hdci

Conclusions:

the estimated intercept (expected presence/absence of Uta lizards when the perimeter to area ratio of an island is 0) is 4.41. This is the median of the posterior distribution for this parameter. When the perimeter to area ratio is 0, Uta lizards are 82.53 times more likely to be present than absent. We can also backtransform to the probability scale. When the perimeter to area ratio is 0, the probability of Uta lizards being present is 0.99
the estimated slope (rate at which presence/absence changes per 1 unit change in perimeter to area ratio), is -0.27 (mean) or (median). The 95% credibility intervals indicate that we are 95% confident that the slope is between -0.52 and -0.08 - e.g. there is a significant positive trend. When backtransformed onto the odds scale, we see that for a one unit increase in perimeter to area ratio, the odds of Uta lizards being present declines by 0.77

bayesplot

polis.rstanarm3 %>% plot(plotfun='mcmc_intervals')

tidy_draws

This is purely a graphical depiction on the posteriors.

polis.rstanarm3 %>% tidy_draws()

spread_draws

polis.rstanarm3 %>% spread_draws(`(Intercept)`, RATIO)

# OR via regex
polis.rstanarm3 %>% spread_draws(`.Intercept.*|RATIO.*`,  regex=TRUE)

posterior_samples

polis.rstanarm3 %>% posterior_samples() %>% as_tibble()

\(R^2\)

polis.rstanarm3 %>% bayes_R2() %>% median_hdci

brms

brms captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

summary(polis.brm3)

##  Family: binomial 
##   Links: mu = logit 
## Formula: PA ~ RATIO 
##    Data: polis (Number of observations: 19) 
## Samples: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup samples = 2400
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     4.57      2.01     1.37     9.10 1.00     2191     2071
## RATIO        -0.28      0.12    -0.55    -0.09 1.00     2084     2214
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Conclusions:

in the Model Info, we are informed that the total MCMC posterior sample size is 2400 and that there were 19 raw observations.
the estimated intercept (expected presence/absence when perimeter to area ratio is 0) is 4.57. This is the mean of the posterior distribution for this parameter. If we back-transform this to the odds-ratio scale, this becomes 96.29. When the perimeter to area ratio is 0, Uta lizards are 96.29 times more likely to be present than absent. We can also backtransform to the probability scale. When the perimeter to area ratio is 0, the probability of Uta lizards being present is 0.99
the estimated slope (rate at which presence/absence changes per 1 unit change in perimeter to area ratio), is -0.28 (mean) or -0.09 (median) with a standard deviation of 0.12. The 90% credibility intervals indicate that we are 90% confident that the slope is between -0.28 and 1 - e.g. there is a significant positive trend. When backtransformed onto the odds scale, we see that for a one unit increase in perimeter to area ratio, the odds of Uta lizards being present declines by 0.76
Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

tidyMCMC

tidyMCMC(polis.brm3$fit, estimate.method='median',  conf.int=TRUE,  conf.method='HPDinterval',  rhat=TRUE, ess=TRUE)

Conclusions:

the estimated intercept (expected presence/absence when perimeter to area ratio is 0) is 4.29. This is the median of the posterior distribution for this parameter. If we back-transform this to the odds-ratio scale, this becomes 73.04. When the perimeter to area ratio is 0, Uta lizards are 73.04 times more likely to be present than absent. We can also backtransform to the probability scale. When the perimeter to area ratio is 0, the probability of Uta lizards being present is 0.99
the estimated slope (rate at which presence/absence changes per 1 unit change in perimeter to area ratio), is -0.26 (median) with a standard error of 0.12. The 95% credibility intervals indicate that we are 95% confident that the slope is between -0.51 and -0.07 - e.g. there is a significant positive trend. When backtransformed onto the odds scale, we see that for a one unit increase in perimeter to area ratio, the odds of Uta lizards being present declines by 0.77
Rhat and number of effective samples for each parameter are also provided and all look good.

gather_draws

polis.draw <- polis.brm3 %>% gather_draws(b_Intercept, b_RATIO)
## OR via regex
polis.draw <- polis.brm3 %>% gather_draws(`b_.*`,  regex=TRUE)
polis.draw

We can then summarise this

polis.draw %>% median_hdci

polis.brm3 %>% 
  gather_draws(b_Intercept, b_RATIO) %>% 
  ggplot() + 
  stat_halfeye(aes(x=.value,  y=.variable)) +
  facet_wrap(~.variable, scales='free')

We could alternatively express the parameters on the odds (odds ratio) scale.

polis.brm3 %>% 
  gather_draws(b_Intercept, b_RATIO) %>%
  group_by(.variable) %>%
  mutate(.value=exp(.value)) %>%
  median_hdci

Conclusions:

the estimated intercept (expected presence/absence of Uta lizards when the perimeter to area ratio of an island is 0) is 4.29. This is the median of the posterior distribution for this parameter. When the perimeter to area ratio is 0, Uta lizards are 73.04 times more likely to be present than absent. We can also backtransform to the probability scale. When the perimeter to area ratio is 0, the probability of Uta lizards being present is 0.99
the estimated slope (rate at which presence/absence changes per 1 unit change in perimeter to area ratio), is -0.26 (mean) or (median). The 95% credibility intervals indicate that we are 95% confident that the slope is between -0.51 and -0.07 - e.g. there is a significant positive trend. When backtransformed onto the odds scale, we see that for a one unit increase in perimeter to area ratio, the odds of Uta lizards being present declines by 0.77

bayesplot

polis.brm3 %>% mcmc_plot(type='intervals')

tidy_draws

This is purely a graphical depiction on the posteriors.

polis.brm3 %>% tidy_draws()

spread_draws

polis.brm3 %>% spread_draws(b_Intercept, b_RATIO)

# OR via regex
polis.brm3 %>% spread_draws(`b_.*`,  regex=TRUE)

posterior_samples

polis.brm3 %>% posterior_samples() %>% as_tibble()

\(R^2\)

polis.brm3 %>% bayes_R2(summary=FALSE) %>% median_hdci

Further analyses

In some disciplines it is useful to be able to calculate an LD50. This is the value along the x-axis that corresponds to a proability of 50% - e.g. the switch-over point in Island perimerter to area Ratio at which the lizards go from more likely to be present to more likely to be absent. It is the inflection point.

It is also the point at which the slope (when backtransformed) is at its steapest and can be calculated as:

\[ -\frac{Intercept}{Slope} \]

rstanarm

polis.rstanarm3 %>% tidy_draws() %>%
  mutate(LD50 = -1*`(Intercept)`/RATIO) %>%
  pull(LD50) %>%
  median_hdci

brms

polis.brm3 %>% tidy_draws() %>%
  mutate(LD50 = -1*b_Intercept/b_RATIO) %>%
  pull(LD50) %>%
  median_hdci

Summary figure

rstanarm

## Using emmeans
polis.grid = with(polis, list(RATIO = seq(min(RATIO), max(RATIO), len=100)))

newdata = emmeans(polis.rstanarm3, ~RATIO, at=polis.grid, type='response') %>% as.data.frame
head(newdata)

ggplot(newdata, aes(y=prob, x=RATIO)) + 
geom_point(data=polis, aes(y=PA)) +
geom_line() + 
geom_ribbon(aes(ymin=lower.HPD, ymax=upper.HPD), fill='blue', alpha=0.3) +
scale_y_continuous('PA') +
scale_x_continuous('RATIO') +
    theme_classic()

brms

## Using emmeans
polis.grid = with(polis, list(RATIO = seq(min(RATIO), max(RATIO), len=100)))

newdata = emmeans(polis.brm3, ~RATIO, at=polis.grid, type='response') %>% as.data.frame
head(newdata)

ggplot(newdata, aes(y=prob, x=RATIO)) + 
geom_point(data=polis, aes(y=PA)) +
geom_line() + 
geom_ribbon(aes(ymin=lower.HPD, ymax=upper.HPD), fill='blue', alpha=0.3) +
scale_y_continuous('PA') +
scale_x_continuous('RATIO') +
    theme_classic()

References

Polis, G. A., S. D. Hurd, C. D. Jackson, and F. Sanchez-Piñero. 1998. “Multifactor Population Limitation: Variable Spatial and Temporal Control of Spiders on Gulf of California Islands.” Ecology 79: 490–502.

Bayesian GLM Part2

Murray Logan

26 December, 2020

Preparations

Scenario

Read in the data

Exploratory data analysis

Fit the model

frequentist

rstanarm

Using default priors

Defining priors

Plotting prior and posterior

brms

Using default priors

Assessing priors

Defining priors

Plotting prior and posterior

Exploring the stan code

MCMC sampling diagnostics

rstanarm

bayesplot

stan plots

ggmcmc

brms

bayesplot

stan plots

ggmcmc

Model validation

rstanarm

pp check

DHARMa residuals

brms

pp check

DHARMa residuals

Partial effects plots

rstanarm

ggpredict

ggemmeans

fitted_draws

brms

conditional_effects

ggpredict

ggemmeans

fitted_draws

Model investigation

rstanarm

summary

tidyMCMC

gather_draws

bayesplot

tidy_draws

spread_draws

posterior_samples

\(R^2\)

brms

summary

tidyMCMC

gather_draws

bayesplot

tidy_draws

spread_draws

posterior_samples

\(R^2\)

Further analyses

rstanarm

brms

Summary figure

rstanarm

brms

References